Experimental enforcement of sex-limited adaptation

knitr::opts_chunk$set(echo = TRUE, warning = FALSE, message = FALSE, cache = TRUE)

Setup for analysis

Load packages

library(tidyverse) # tidy style coding
library(brms) # Bayesian models
library(loo) # to use information criteria in brms models
library(tidybayes) # Bayesian aesthetics
library(MetBrewer) # colours
library(kableExtra) # tables
library(patchwork) # putting plots together
library(DT) # for search- and saveable tables
library(pander) # for simpler tables

Load in the data

\(~\)

fitness_data <- read_csv("data/SLC_fitness_data.csv") %>% 
  mutate(Fitness_vial_ID = as.factor(Fitness_vial_ID),
         Block = as.factor(Block),
         Population = as.factor(Population),
         Treatment = as.factor(Treatment),
         GFP = as.factor(GFP),
         Sex = as.factor(Sex),
         Rearing_vial = as.factor(Rearing_vial),
         Total_red_offspring = Red_female_offspring + Red_male_offspring,
         Total_bw_offspring = bw_female_offspring + bw_male_offspring,
         Total_offspring = Total_red_offspring + Total_bw_offspring) %>% 
  rename(Evolutionary_treatment = Treatment)

# Create a function to build HTML searchable tables

my_data_table <- function(df){
  datatable(
    df, rownames=FALSE,
    autoHideNavigation = TRUE,
    extensions = c("Scroller",  "Buttons"),
    options = list(
      dom = 'Bfrtip',
      deferRender=TRUE,
      scrollX=TRUE, scrollY=400,
      scrollCollapse=TRUE,
      buttons =
        list('pageLength', 'colvis', 'csv', list(
          extend = 'pdf',
          pageSize = 'A4',
          orientation = 'landscape',
          filename = 'fitness_data')),
      pageLength = 692
    )
  )
}


my_data_table(fitness_data %>% select(-Comment))

Column explanations

Fitness_vial_ID: a unique identifier for each trial of the fitness assay.

Block: the experiment was run in three distinct blocks, using flies from separate generations.

Population: we measured the fitness of flies from 12 independent populations that contained autosomes that had undergone experimental evolution.

Evolutionary_treatment: the populations carried autosomes that had been exposed to one of three evolutionary conditions for 20 generations: a female-limited response to selection, a male-limited response and a control condition where the evolutionary response was unconstrained.

GFP: the GFP marker carried by the population. UBI indicates the presence of a transgene that encodes ubiquitous expression of GFP, while 3xP indicates the presence of a different transgene that encodes the expression of GFP in the ocelli.

Sex: the sex of the individuals that we were measuring the fitness of.

Rearing_vial: the vial the treatment flies used in the trial developed in. This variable is included to capture variation explained by the rearing environment e.g. small differences in food moisture content or quantity. Note that females and males can have the same rearing vial as the sexes were reared together.

Red_female_offspring: the number of adult female offspring sired/produced by flies sourced from one of the 12 populations.

Red_male_offspring: the number of adult male offspring sired/produced by flies sourced from one of the 12 populations.

bw_female_offspring: the number of adult female offspring sired/produced by the competitor flies in our fitness assay. bw is a recessive allele that encodes brown eye colour.

bw_male_offspring: the number of adult male offspring sired/produced by the competitor flies in our fitness assay. bw is a recessive allele that encodes brown eye colour.

Total_red_offspring: the total number (sexes pooled) of adult offspring sired/produced by flies sourced from one of the 12 populations.

Total_bw_offspring: the total number (sexes pooled) of adult offspring sired/produced by competitor flies.

Total_offspring: the total number (sexes and eye colours pooled) of adult offspring counted in each vial.

\(~\)

Modelling approach

\(~\)

Female and male fitness are fundamentally different concepts / traits. There are also several differences between our female and male fitness assays. The major difference is that the male assay contains half the number of females in any given vial than does the female assay. The logic behind this design choice is that sexually selected processes are a more important determinant of male fitness than they are female fitness, so any fitness differences may only be observed when competition for fertilisations is high.

For these reasons, we choose to split the data up and model female and male fitness separately.

female_fitness <- 
  fitness_data %>%
  filter(Sex == "Female")

male_fitness <- 
  fitness_data %>%
  filter(Sex == "Male") %>% 
  mutate(prop_red = Total_red_offspring / Total_offspring)

We fit the following fixed and random effects to model female and male fitness. Our aim is to estimate the causal effect that evolutionary treatment has on fitness.

Fixed effects

Block: fitness might differ between the three distinct blocks we split our experiment up into. Blocks differed temporally, used flies from different generations and different batches of food. It is also possible that there were minor fluctuations in the lighting and temperature environment experienced during development between blocks. Each of these variables may introduce variation into our fitness measurements, that can be accounted for by including the Block variable in our model.

GFP: it is possible that fitness may be affected by the GFP transgene carried by each population. For example, one could imagine that any unintended fitness effects of a transgene might be of greater magnitude if it is expressed in a larger proportion of tissues, as is the case for the UBI transgene versus the 3xP transgene. Note that each GFP type is carried by an equal number of populations from each of the three evolutionary treatments.

Evolutionary_treatment: this is the inheritance regime that the autosomes carried by each of the populations were subject to for 20 generations. There are three levels: populations carrying female-adapted autosomes, populations containing male-adapted autosomes and populations carrying control autosomes that experienced an unmanipulated inheritance regime. We are designing our model to test for a causal effect of this variable.

Random effects

Rearing_vial: the vial individual flies developed within may introduce further variation into our response variable. Like Block this variable controls for micro-environmental variation.

Population: our design contained 12 independent populations of autosomes that originated from a single outbred laboratory fly population. The populations were split and autosomes from each were subjected to one of the three evolution treatments for 20 generations. 4 populations experienced a female-limited inheritance regime, 4 a male-limited regime and 4 an unlimited or control regime.

\(~\)

Female fitness

\(~\)

To estimate the fitness of females carrying each of the three autosome types, we placed three experimental females into a yeasted vial with three female competitors that carried the bw mutation. We then introduced six males that also carried the bw mutation. We allowed them to mate and oviposit for three days, then removed all adults from the vial. 12 days later we counted all of the adult progeny in the vial and scored them for eye-colour. Progeny with red eyes were produced by the experimental females ( bw is recessive) and progeny with brown eyes were produced by the competitor females. We calculated fitness as the proportion of red eyed offspring in the vial.

\(~\)

Modelling

female_fitness.01 <-
  brm(Total_red_offspring | trials(Total_offspring) ~ 1 + Evolutionary_treatment + GFP + Block + (1|Population) + (1|Rearing_vial),
      data = female_fitness, family = binomial(), 
      prior = c(prior(normal(0, 1.5), class = Intercept),
                prior(normal(0, 2), class = b),
                prior(exponential(1), class = sd)),
      iter = 10000, warmup = 8000, chains = 4, cores = 4,
      control = list(adapt_delta = 0.999, max_treedepth = 15),
      seed = 2, file = "Fits/female_fitness.01")

print(female_fitness.01)

##  Family: binomial 
##   Links: mu = logit 
## Formula: Total_red_offspring | trials(Total_offspring) ~ 1 + Evolutionary_treatment + GFP + Block + (1 | Population) + (1 | Rearing_vial) 
##    Data: female_fitness (Number of observations: 327) 
##   Draws: 4 chains, each with iter = 10000; warmup = 8000; thin = 1;
##          total post-warmup draws = 8000
## 
## Group-Level Effects: 
## ~Population (Number of levels: 12) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     0.06      0.05     0.00     0.18 1.01     1090     2583
## 
## ~Rearing_vial (Number of levels: 70) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     0.27      0.03     0.22     0.33 1.00     2030     3600
## 
## Population-Level Effects: 
##                                      Estimate Est.Error l-95% CI u-95% CI Rhat
## Intercept                                1.09      0.09     0.91     1.27 1.00
## Evolutionary_treatmentFemale_limited     0.20      0.10    -0.00     0.40 1.00
## Evolutionary_treatmentMale_limited       0.16      0.10    -0.02     0.36 1.00
## GFPUBI                                  -0.32      0.08    -0.48    -0.16 1.00
## Block2                                  -0.58      0.08    -0.75    -0.42 1.00
## Block3                                  -0.61      0.08    -0.77    -0.45 1.00
##                                      Bulk_ESS Tail_ESS
## Intercept                                2314     3772
## Evolutionary_treatmentFemale_limited     2500     3041
## Evolutionary_treatmentMale_limited       2884     3849
## GFPUBI                                   2625     3724
## Block2                                   2042     2928
## Block3                                   1996     2892
## 
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).

Lets see if we have any points that have a strong influence on the posterior distribution

#female_fitness.01 <- add_criterion(female_fitness.01, criterion = "loo", file = "fits/female_fitness.01")

#female_fitness.01 <- add_criterion(female_fitness.01, criterion = "waic", file = "fits/female_fitness.01")

loo(female_fitness.01)

## 
## Computed from 8000 by 327 log-likelihood matrix
## 
##          Estimate    SE
## elpd_loo  -1913.3  73.7
## p_loo       302.7  23.3
## looic      3826.6 147.4
## ------
## Monte Carlo SE of elpd_loo is NA.
## 
## Pareto k diagnostic values:
##                          Count Pct.    Min. n_eff
## (-Inf, 0.5]   (good)     241   73.7%   565       
##  (0.5, 0.7]   (ok)        53   16.2%   138       
##    (0.7, 1]   (bad)       26    8.0%   13        
##    (1, Inf)   (very bad)   7    2.1%   3         
## See help('pareto-k-diagnostic') for details.

There appear to be many points that are highly influential on the posterior. Why? It appears that this is explained by small sample sizes for some of the rearing vials. For example, female flies reared in vial 24 are only used in a single fitness assay, making estimation for this vial very difficult. It also means that the single observation for that vial has a huge effect on its posterior distribution, which returns a very high pareto k value.

There are several ways to limit the influence of these ‘outliers’ on the posterior. The first is to remove Rearing_vial from the model. While not run here, this reduces the number of highly influential points to 0. Alternatively, we can specify a different distribution family that is better equipped to deal with extreme values at the tails. This allows us to account for the potentially important variation coded by Rearing_vial. A more flexible version of the binomial is the betabinomial.

However, the beta-binomial is not a native family in brms, we need to create the distribution family using the custom_family() function. The code below is taken directly from the brms_customfamilies vignette, which can be viewed here.

beta_binomial2 <- custom_family(
  "beta_binomial2", dpars = c("mu", "phi"),
  links = c("logit", "log"), lb = c(NA, 0),
  type = "int", vars = "vint1[n]"
)


stan_funs <- "
  real beta_binomial2_lpmf(int y, real mu, real phi, int T) {
    return beta_binomial_lpmf(y | T, mu * phi, (1 - mu) * phi);
  }
  int beta_binomial2_rng(real mu, real phi, int T) {
    return beta_binomial_rng(T, mu * phi, (1 - mu) * phi);
  }
"

stanvars <- stanvar(scode = stan_funs, block = "functions")

Below we fit the same model, but this time specify a beta-binomial distribution.

female_fitness.02 <-
  brm(Total_red_offspring | vint(Total_offspring) ~ 1 + Evolutionary_treatment + GFP + Block + (1|Population) + (1|Rearing_vial),
      data = female_fitness, family = beta_binomial2, 
      prior = c(prior(normal(0, 1.5), class = Intercept),
                prior(normal(0, 2), class = b),
                prior(exponential(1), class = sd)),
      iter = 10000, warmup = 8000, chains = 4, cores = 4,
      control = list(adapt_delta = 0.95, max_treedepth = 10),
      seed = 2, stanvars = stanvars, file = "Fits/female_fitness.02")

print(female_fitness.02)

##  Family: beta_binomial2 
##   Links: mu = logit; phi = identity 
## Formula: Total_red_offspring | vint(Total_offspring) ~ 1 + Evolutionary_treatment + GFP + Block + (1 | Population) + (1 | Rearing_vial) 
##    Data: female_fitness (Number of observations: 327) 
##   Draws: 4 chains, each with iter = 10000; warmup = 8000; thin = 1;
##          total post-warmup draws = 8000
## 
## Group-Level Effects: 
## ~Population (Number of levels: 12) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     0.08      0.05     0.00     0.20 1.00     1658     2882
## 
## ~Rearing_vial (Number of levels: 70) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     0.17      0.04     0.09     0.25 1.00     2319     3274
## 
## Population-Level Effects: 
##                                      Estimate Est.Error l-95% CI u-95% CI Rhat
## Intercept                                1.13      0.09     0.94     1.31 1.00
## Evolutionary_treatmentFemale_limited     0.17      0.10    -0.02     0.37 1.00
## Evolutionary_treatmentMale_limited       0.18      0.10    -0.01     0.38 1.00
## GFPUBI                                  -0.33      0.08    -0.49    -0.16 1.00
## Block2                                  -0.64      0.08    -0.79    -0.48 1.00
## Block3                                  -0.63      0.07    -0.78    -0.49 1.00
##                                      Bulk_ESS Tail_ESS
## Intercept                                4700     5092
## Evolutionary_treatmentFemale_limited     4537     5200
## Evolutionary_treatmentMale_limited       4465     4696
## GFPUBI                                   4853     5113
## Block2                                   5339     5608
## Block3                                   4865     5524
## 
## Family Specific Parameters: 
##     Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## phi    33.64      3.60    27.06    41.11 1.00     6092     5523
## 
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).

#female_fitness.02 <- add_criterion(female_fitness.02, criterion = "loo", file = "fits/female_fitness.02")

#female_fitness.02 <- add_criterion(female_fitness.02, criterion = "waic", file = "fits/female_fitness.02")

We need to write some additional code to get some post processing stuff i.e. LOO to work. Code courtesy of the brms_customfamilies vignette, which can be viewed here.

# we need to write some additional code to get some post-processing stuff to work

expose_functions(female_fitness.02, vectorize = TRUE)


log_lik_beta_binomial2 <- function(i, prep) {
  mu <- brms::get_dpar(prep, "mu", i = i)
  phi <- brms::get_dpar(prep, "phi", i = i)
  trials <- prep$data$vint1[i]
  y <- prep$data$Y[i]
  beta_binomial2_lpmf(y, mu, phi, trials)
}

posterior_predict_beta_binomial2 <- function(i, prep, ...) {
  mu <- brms::get_dpar(prep, "mu", i = i)
  phi <- brms::get_dpar(prep, "phi", i = i)
  trials <- prep$data$vint1[i]
  beta_binomial2_rng(mu, phi, trials)
}

posterior_epred_beta_binomial2 <- function(prep) {
  mu <- brms::get_dpar(prep, "mu")
  trials <- prep$data$vint1
  trials <- matrix(trials, nrow = nrow(mu), ncol = ncol(mu), byrow = TRUE)
  mu * trials
}

Run LOO to see if we’ve handled the influential points and then rank models with the loo_compare() function.

loo(female_fitness.02)

## 
## Computed from 8000 by 327 log-likelihood matrix
## 
##          Estimate   SE
## elpd_loo  -1395.5 13.0
## p_loo        36.8  3.3
## looic      2791.0 26.0
## ------
## Monte Carlo SE of elpd_loo is 0.1.
## 
## Pareto k diagnostic values:
##                          Count Pct.    Min. n_eff
## (-Inf, 0.5]   (good)     326   99.7%   1299      
##  (0.5, 0.7]   (ok)         1    0.3%   521       
##    (0.7, 1]   (bad)        0    0.0%   <NA>      
##    (1, Inf)   (very bad)   0    0.0%   <NA>      
## 
## All Pareto k estimates are ok (k < 0.7).
## See help('pareto-k-diagnostic') for details.

loo_compare(female_fitness.01, female_fitness.02, criterion = "waic")

##                   elpd_diff se_diff
## female_fitness.02    0.0       0.0 
## female_fitness.01 -505.8      62.7

The beta-binomial model looks good. It fixes our influential points problem and is preferred by the LOO metric.

Conduct a posterior predictive check to confirm our model is doing what we want it to.

pp_check(female_fitness.02, type = "hist", ndraws = 11, binwidth = 10) +
  theme_minimal() +
  theme(panel.background = element_blank())

The posterior predictive distribution matches our raw data quite well, indicating the model is functioning as we wanted.

\(~\)

Derive predictions from the posterior

Get predictions from the model and present them in a Table

draws <- as_draws_df(female_fitness.02)

draws_female <-
  draws  %>% 
  mutate(`Female-limited` = inv_logit_scaled(b_Intercept + b_Evolutionary_treatmentFemale_limited),
         `Male-limited` = inv_logit_scaled(b_Intercept + b_Evolutionary_treatmentMale_limited),
         Control = inv_logit_scaled(b_Intercept)) %>% 
  select(Control, `Female-limited`, `Male-limited`) %>% 
    pivot_longer(cols = c(Control, `Female-limited`, `Male-limited`),
                 names_to = "Evolutionary_treatment") %>% 
  rename(prop_focal_offspring = value) %>% 
  arrange(Evolutionary_treatment)

draws_female %>% 
  group_by(Evolutionary_treatment) %>% 
  summarise(`Estimated prop. of offspring produced` = median(prop_focal_offspring),
            `2.5%` = quantile(prop_focal_offspring, probs = 0.025),
            `97.5%` = quantile(prop_focal_offspring, probs = 0.975)) %>% 
  rename(`Inheritance treatment` = Evolutionary_treatment) %>% 
  pander(split.cell = 40, split.table = Inf, round = 3)

Inheritance treatment	Estimated prop. of offspring produced	2.5%	97.5%
Control	0.755	0.719	0.787
Female-limited	0.786	0.751	0.816
Male-limited	0.788	0.754	0.817

# now find the differences between the control and the sex-limited treatments

# the inv_logit_scaled() function converts the posterior draws onto the response scale 

  draws %>% 
  mutate(p_control =  inv_logit_scaled(b_Intercept),
         p_female = inv_logit_scaled(b_Evolutionary_treatmentFemale_limited + b_Intercept),
         p_male = inv_logit_scaled(b_Evolutionary_treatmentMale_limited + b_Intercept),
         `Female-limited` = p_female / p_control,
         `Male-limited` = p_male / p_control) %>% 
  gather(key = `difference comparison`, value = `% difference`) %>% 
  filter(`difference comparison` == c("Female-limited", "Male-limited")) %>% 
  group_by(`difference comparison`)  %>% 
  summarise(`Mean proportion of offspring produced relative to control`  = mean(`% difference`),
            `2.5%` = quantile(`% difference`, probs = 0.025),
            `97.5%` = quantile(`% difference`, probs = 0.975)) %>% 
  rename(`Inheritance treatment` = `difference comparison`) %>% 
  pander(split.cell = 40, split.table = Inf, round = 3)

Inheritance treatment	Mean proportion of offspring produced relative to control	2.5%	97.5%
Female-limited	1.041	0.995	1.089
Male-limited	1.043	0.997	1.092

\(~\)

Male fitness

\(~\)

To estimate the fitness of males carrying each of the three chromosome types, we conducted an experiment very similar to the female fitness assay. However, because male fitness has stronger covariance with fertilisation events than does female fitness, we conducted the male fitness assay with a 1:2 sex ratio (female:male) rather than the 1:1 ratio used in the female assay. This increases the strength of sexual selection and is a more appopriate way to expose differences in fitness between groups of males. As with the females, we calculated fitness as the proportion of red eyed offspring in the vial.

Alternative modelling approaches

Binomial model with rearing vial

male_fitness.01 <-
  brm(Total_red_offspring | trials(Total_offspring) ~ 1 + Evolutionary_treatment + GFP + Block + (1|Population) + (1|Rearing_vial),
      data = male_fitness, family = binomial(), 
      prior = c(prior(normal(0, 1.5), class = Intercept),
                prior(normal(0, 2), class = b),
                prior(exponential(1), class = sd)),
      iter = 10000, warmup = 8000, chains = 4, cores = 4,
      control = list(adapt_delta = 0.999, max_treedepth = 15),
      seed = 2, file = "Fits/male_fitness.01")


print(male_fitness.01)

##  Family: binomial 
##   Links: mu = logit 
## Formula: Total_red_offspring | trials(Total_offspring) ~ 1 + Evolutionary_treatment + GFP + Block + (1 | Population) + (1 | Rearing_vial) 
##    Data: male_fitness (Number of observations: 360) 
##   Draws: 4 chains, each with iter = 10000; warmup = 8000; thin = 1;
##          total post-warmup draws = 8000
## 
## Group-Level Effects: 
## ~Population (Number of levels: 12) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     0.30      0.15     0.06     0.64 1.00     1307     2009
## 
## ~Rearing_vial (Number of levels: 72) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     0.61      0.06     0.50     0.75 1.00     2567     4428
## 
## Population-Level Effects: 
##                                      Estimate Est.Error l-95% CI u-95% CI Rhat
## Intercept                                0.66      0.26     0.13     1.18 1.00
## Evolutionary_treatmentFemale_limited     0.13      0.29    -0.44     0.73 1.00
## Evolutionary_treatmentMale_limited       0.27      0.30    -0.34     0.86 1.00
## GFPUBI                                  -0.08      0.23    -0.57     0.35 1.00
## Block2                                   0.92      0.18     0.56     1.28 1.00
## Block3                                   0.05      0.18    -0.30     0.40 1.00
##                                      Bulk_ESS Tail_ESS
## Intercept                                5006     4605
## Evolutionary_treatmentFemale_limited     5510     4532
## Evolutionary_treatmentMale_limited       5350     4628
## GFPUBI                                   5054     5100
## Block2                                   3643     4695
## Block3                                   3626     4385
## 
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).

Like with the female fitness model, lets see if we have any points that have a strong influence on the posterior distribution. Run LOO:

#male_fitness.01 <- add_criterion(male_fitness.04, criterion = "loo", file = "fits/male_fitness.04")
#male_fitness.01 <- add_criterion(male_fitness.01, criterion = "waic", file = "fits/male_fitness.01")

loo(male_fitness.01)

## 
## Computed from 8000 by 360 log-likelihood matrix
## 
##          Estimate    SE
## elpd_loo  -5494.1 310.1
## p_loo      1393.9 101.3
## looic     10988.2 620.3
## ------
## Monte Carlo SE of elpd_loo is NA.
## 
## Pareto k diagnostic values:
##                          Count Pct.    Min. n_eff
## (-Inf, 0.5]   (good)     141   39.2%   218       
##  (0.5, 0.7]   (ok)        75   20.8%   66        
##    (0.7, 1]   (bad)       65   18.1%   15        
##    (1, Inf)   (very bad)  79   21.9%   1         
## See help('pareto-k-diagnostic') for details.

Conduct the posterior predictive check…

pp_check(male_fitness.01, type = "hist", ndraws = 11, binwidth = 10) +
  theme_minimal() +
  theme(panel.background = element_blank())

Binomial model without rearing vial

Like with the female model, there appears to be many points that are highly influential on the posterior. I remove Rearing_vial once more to see if this is the variable responsible for making certain points highly influential.

male_fitness.02 <-
  brm(Total_red_offspring | trials(Total_offspring) ~ 1 + Evolutionary_treatment + GFP + Block + (1|Population),
      data = male_fitness, family = binomial(), 
      prior = c(prior(normal(0, 1.5), class = Intercept),
                prior(normal(0, 2), class = b),
                prior(exponential(1), class = sd)),
      iter = 4000, warmup = 2000, chains = 4, cores = 4,
      control = list(adapt_delta = 0.8, max_treedepth = 10),
      seed = 2, file = "Fits/male_fitness.02")

Run LOO

#male_fitness.02 <- add_criterion(male_fitness.02, criterion = "loo", file = "fits/male_fitness.02")

#male_fitness.02 <- add_criterion(male_fitness.02, criterion = "waic", file = "fits/male_fitness.02")

loo(male_fitness.02)

## 
## Computed from 8000 by 360 log-likelihood matrix
## 
##          Estimate    SE
## elpd_loo  -5624.8 304.2
## p_loo       347.0  27.0
## looic     11249.5 608.4
## ------
## Monte Carlo SE of elpd_loo is NA.
## 
## Pareto k diagnostic values:
##                          Count Pct.    Min. n_eff
## (-Inf, 0.5]   (good)     316   87.8%   372       
##  (0.5, 0.7]   (ok)        27    7.5%   124       
##    (0.7, 1]   (bad)       15    4.2%   25        
##    (1, Inf)   (very bad)   2    0.6%   13        
## See help('pareto-k-diagnostic') for details.

Conduct the posterior predictive check…

pp_check(male_fitness.02, type = "hist", ndraws = 11, binwidth = 10) +
  theme_minimal() +
  theme(panel.background = element_blank())

Removing Rearing_vial reduces the number of highly influential points, however this time 14 remain…

Beta-binomial model

male_fitness.03 <-
  brm(Total_red_offspring | vint(Total_offspring) ~ 1 + Evolutionary_treatment + GFP + Block + (1|Population) + (1|Rearing_vial),
      data = male_fitness, family = beta_binomial2, 
      prior = c(prior(normal(0, 1.5), class = Intercept),
                prior(normal(0, 2), class = b),
                prior(exponential(1), class = sd)),
      iter = 10000, warmup = 8000, chains = 4, cores = 4,
      stanvars = stanvars,
      control = list(adapt_delta = 0.9, max_treedepth = 10),
      seed = 2, file = "Fits/male_fitness.03")

print(male_fitness.03)

##  Family: beta_binomial2 
##   Links: mu = logit; phi = identity 
## Formula: Total_red_offspring | vint(Total_offspring) ~ 1 + Evolutionary_treatment + GFP + Block + (1 | Population) + (1 | Rearing_vial) 
##    Data: male_fitness (Number of observations: 360) 
##   Draws: 4 chains, each with iter = 10000; warmup = 8000; thin = 1;
##          total post-warmup draws = 8000
## 
## Group-Level Effects: 
## ~Population (Number of levels: 12) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     0.29      0.13     0.06     0.60 1.00     2214     2085
## 
## ~Rearing_vial (Number of levels: 72) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     0.16      0.10     0.01     0.38 1.00     2143     3826
## 
## Population-Level Effects: 
##                                      Estimate Est.Error l-95% CI u-95% CI Rhat
## Intercept                                0.77      0.23     0.31     1.22 1.00
## Evolutionary_treatmentFemale_limited     0.42      0.27    -0.11     0.95 1.00
## Evolutionary_treatmentMale_limited       0.10      0.26    -0.42     0.63 1.00
## GFPUBI                                  -0.34      0.22    -0.78     0.10 1.00
## Block2                                   0.98      0.15     0.69     1.28 1.00
## Block3                                  -0.03      0.14    -0.30     0.25 1.00
##                                      Bulk_ESS Tail_ESS
## Intercept                                5593     5203
## Evolutionary_treatmentFemale_limited     5048     4596
## Evolutionary_treatmentMale_limited       5141     5227
## GFPUBI                                   5347     4793
## Block2                                   9012     5026
## Block3                                   9113     5995
## 
## Family Specific Parameters: 
##     Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## phi     3.65      0.31     3.08     4.29 1.00    10770     5671
## 
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).

Quick bit of code to get our custom family post-processing to work

expose_functions(male_fitness.03, vectorize = TRUE)

Run LOO

#male_fitness.03 <- add_criterion(male_fitness.03, criterion = "loo", file = "fits/male_fitness.03")

#male_fitness.03 <- add_criterion(male_fitness.03, criterion = "waic", file = "fits/male_fitness.03")

loo(male_fitness.03)

## 
## Computed from 8000 by 360 log-likelihood matrix
## 
##          Estimate   SE
## elpd_loo  -1561.9 21.4
## p_loo        23.2  1.9
## looic      3123.7 42.9
## ------
## Monte Carlo SE of elpd_loo is NA.
## 
## Pareto k diagnostic values:
##                          Count Pct.    Min. n_eff
## (-Inf, 0.5]   (good)     358   99.4%   2421      
##  (0.5, 0.7]   (ok)         1    0.3%   2668      
##    (0.7, 1]   (bad)        0    0.0%   <NA>      
##    (1, Inf)   (very bad)   1    0.3%   4000      
## See help('pareto-k-diagnostic') for details.

loo_compare(male_fitness.01, male_fitness.02, male_fitness.03, criterion = "loo")

##                 elpd_diff se_diff
## male_fitness.03     0.0       0.0
## male_fitness.01 -3932.2     305.5
## male_fitness.02 -4062.9     298.9

There is only a single point now that LOO is worried about, even while including Rearing_vial as a random effect. The beta-binomial model also performs the best by the LOO and WAIC metrics.

Conduct the posterior predictive check…

pp_check(male_fitness.03, type = "hist", ndraws = 11, binwidth = 10) +
  theme_minimal() +
  theme(panel.background = element_blank())

Zero-one-inflated beta

Another interesting aspect to these data are the prevalence of 1 values. That is, there are many trials where the red-eyed males carrying the experimentally evolved autosomes sire all of the offspring produced in the vial. These extreme values are not handled well by binomial models. Instead we can model the proportion of red-eyed offspring using a zero-one-inflated beta model.

male_fitness.04 <-
  brm(prop_red ~ 1 + Evolutionary_treatment + GFP + Block + (1|Population) + (1|Rearing_vial),
      data = male_fitness, family = zero_one_inflated_beta(), 
      prior = c(prior(normal(0, 1.5), class = Intercept),
                prior(normal(0, 2), class = b),
                prior(exponential(1), class = sd)),
      iter = 10000, warmup = 8000, chains = 4, cores = 4,
      control = list(adapt_delta = 0.999, max_treedepth = 15),
      seed = 2, file = "Fits/male_fitness.04")


print(male_fitness.01)

##  Family: binomial 
##   Links: mu = logit 
## Formula: Total_red_offspring | trials(Total_offspring) ~ 1 + Evolutionary_treatment + GFP + Block + (1 | Population) + (1 | Rearing_vial) 
##    Data: male_fitness (Number of observations: 360) 
##   Draws: 4 chains, each with iter = 10000; warmup = 8000; thin = 1;
##          total post-warmup draws = 8000
## 
## Group-Level Effects: 
## ~Population (Number of levels: 12) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     0.30      0.15     0.06     0.64 1.00     1307     2009
## 
## ~Rearing_vial (Number of levels: 72) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     0.61      0.06     0.50     0.75 1.00     2567     4428
## 
## Population-Level Effects: 
##                                      Estimate Est.Error l-95% CI u-95% CI Rhat
## Intercept                                0.66      0.26     0.13     1.18 1.00
## Evolutionary_treatmentFemale_limited     0.13      0.29    -0.44     0.73 1.00
## Evolutionary_treatmentMale_limited       0.27      0.30    -0.34     0.86 1.00
## GFPUBI                                  -0.08      0.23    -0.57     0.35 1.00
## Block2                                   0.92      0.18     0.56     1.28 1.00
## Block3                                   0.05      0.18    -0.30     0.40 1.00
##                                      Bulk_ESS Tail_ESS
## Intercept                                5006     4605
## Evolutionary_treatmentFemale_limited     5510     4532
## Evolutionary_treatmentMale_limited       5350     4628
## GFPUBI                                   5054     5100
## Block2                                   3643     4695
## Block3                                   3626     4385
## 
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).

#male_fitness.04 <- add_criterion(male_fitness.04, criterion = "loo", file = "fits/male_fitness.04")
#male_fitness.01 <- add_criterion(male_fitness.01, criterion = "waic", file = "fits/male_fitness.01")

loo(male_fitness.04)

## 
## Computed from 8000 by 359 log-likelihood matrix
## 
##          Estimate   SE
## elpd_loo    -28.7 17.9
## p_loo        16.7  1.4
## looic        57.4 35.9
## ------
## Monte Carlo SE of elpd_loo is 0.0.
## 
## All Pareto k estimates are good (k < 0.5).
## See help('pareto-k-diagnostic') for details.

Note that we no longer have any values with high influence.

Conduct the posterior predictive check…

pp_check(male_fitness.04, type = "hist", ndraws = 11, binwidth = 0.1) +
  theme_minimal() +
  theme(panel.background = element_blank())

Derive estimates from posterior

Get predictions from the model and present them in a Table

draws_m <- as_draws_df(male_fitness.03)

draws_male <-
  draws_m  %>% 
  mutate(`Female-limited` = inv_logit_scaled(b_Intercept + b_Evolutionary_treatmentFemale_limited),
         `Male-limited` = inv_logit_scaled(b_Intercept + b_Evolutionary_treatmentMale_limited),
         Control = inv_logit_scaled(b_Intercept)) %>% 
  select(Control, `Female-limited`, `Male-limited`) %>% 
    pivot_longer(cols = c(Control, `Female-limited`, `Male-limited`),
                 names_to = "Evolutionary_treatment") %>% 
  rename(prop_focal_offspring = value)

draws_male %>% 
  group_by(Evolutionary_treatment) %>% 
  summarise(`Estimated prop. of offspring sired` = median(prop_focal_offspring),
            `2.5%` = quantile(prop_focal_offspring, probs = 0.025),
            `97.5%` = quantile(prop_focal_offspring, probs = 0.975)) %>% 
  rename(`Inheritance treatment` = Evolutionary_treatment) %>% 
  pander(split.cell = 40, split.table = Inf, round = 3)

Inheritance treatment	Estimated prop. of offspring sired	2.5%	97.5%
Control	0.684	0.577	0.772
Female-limited	0.768	0.674	0.838
Male-limited	0.707	0.604	0.792

# now find the differences between the control and the sex-limited Evolution_treatments

# the inv_logit_scaled() function converts the posterior draws onto the response scale 

draws_m %>% 
  mutate(p_control =  inv_logit_scaled(b_Intercept),
         p_female = inv_logit_scaled(b_Evolutionary_treatmentFemale_limited + b_Intercept),
         p_male = inv_logit_scaled(b_Evolutionary_treatmentMale_limited + b_Intercept),
         `Female-limited` = p_female / p_control,
         `Male-limited` = p_male / p_control) %>% 
  gather(key = `difference comparison`, value = `% difference`) %>% 
  filter(`difference comparison` == c("Female-limited", "Male-limited")) %>% 
  group_by(`difference comparison`)  %>% 
  summarise(`Mean proportion of offspring sired relative to control`  = mean(`% difference`),
            `2.5%` = quantile(`% difference`, probs = 0.025),
            `97.5%` = quantile(`% difference`, probs = 0.975)) %>% 
  rename(`Inheritance treatment` = `difference comparison`) %>% 
  pander(split.cell = 40, split.table = Inf, round = 3)

Inheritance treatment	Mean proportion of offspring sired relative to control	2.5%	97.5%
Female-limited	1.126	0.966	1.323
Male-limited	1.035	0.878	1.214

Building Figure 1

# female plots

f_1 <-
  draws_female %>% 
  ggplot(aes(Evolutionary_treatment, prop_focal_offspring)) + 
  stat_eye(aes(fill = Evolutionary_treatment), .width = c(0.66, 0.95), alpha = 1) + # width indicates the uncertainty intervals: here we have 66% and 95% intervals
  scale_fill_manual(values = met.brewer("Hiroshige", 3)) +
  scale_y_continuous("Female fitness\n(proportion of offspring produced)",
                     breaks = c(0.7, 0.8)) +
  #coord_cartesian(ylim = c(0.4, 1)) +
  xlab(NULL) +
  theme_bw() + 
  theme(legend.position = "none",
        panel.grid.minor = element_blank())

f_2 <-
  draws %>% 
  mutate(`Female-limited` = b_Evolutionary_treatmentFemale_limited,
         `Male-limited` = b_Evolutionary_treatmentMale_limited) %>% 
  gather(key = parameter, value = logodds) %>% 
  filter(parameter == c("Female-limited", "Male-limited")) %>%
  as_tibble() %>% 
  
  ggplot(aes(parameter, logodds)) + 
  stat_halfeye(aes(fill = parameter), .width = c(0.66, 0.95)) + # width indicates the uncertainty intervals: here we have 66% and 95% intervals
  scale_fill_manual(values = c("Female-limited" = "#ffd06f", "Male-limited" = "#72bcd5")) + 
  coord_flip() +
  geom_hline(yintercept = 0, linetype = 2) +
  scale_y_continuous(breaks = c(-1, 0, 1)) +
  xlab(NULL) +
  ylab(NULL) +
  theme_bw() + 
  theme(legend.position = "none",
        panel.grid.minor = element_blank())

# male plots

f_3 <-
  draws_male %>% 
  ggplot(aes(Evolutionary_treatment, prop_focal_offspring)) + 
  stat_eye(aes(fill = Evolutionary_treatment), .width = c(0.66, 0.95), alpha = 1) + # width indicates the uncertainty intervals: here we have 66% and 95% intervals
  scale_fill_manual(values = met.brewer("Hiroshige", 3)) +
  scale_y_continuous("Male fitness\n(proportion of offspring sired)",
                     breaks = c(0.4, 0.6, 0.8)) +
  #coord_cartesian(ylim = c(0.4, 1)) +
  xlab("Inheritance treatment") +
  theme_bw() + 
  theme(legend.position = "none",
        panel.grid.minor = element_blank())

f_4 <-
  draws_m %>% 
  mutate(`Female-limited` = b_Evolutionary_treatmentFemale_limited,
         `Male-limited` = b_Evolutionary_treatmentMale_limited) %>% 
  gather(key = parameter, value = logodds) %>% 
  filter(parameter == c("Female-limited", "Male-limited")) %>%
  as_tibble() %>% 
  #mutate(parameter =factor(parameter, levels=c("SD-Mad", "SD-72", "SD-5"))) %>%
  
  ggplot(aes(parameter, logodds)) + 
  stat_halfeye(aes(fill = parameter), .width = c(0.66, 0.95)) + # width indicates the uncertainty intervals: here we have 66% and 95% intervals
  scale_fill_manual(values = c("Female-limited" = "#ffd06f", "Male-limited" = "#72bcd5")) + 
  coord_flip(ylim = c(-1, 1)) +
  geom_hline(yintercept = 0, linetype = 2) +
  scale_y_continuous(breaks = c(-1, 0, 1)) +
  xlab(NULL) +
  ylab("Log-odds mean difference from control") +
  theme_bw() + 
  theme(legend.position = "none",
        panel.grid.minor = element_blank())


(f_1 + f_2) /(f_3 + f_4)

Figure 1: a shows the estimated mean female fitness for flies carrying autosomes that had previously experienced unconstrained inheritance (control), female-limited inheritance or male-limited inheritance. b shows the log-odds mean difference in means between the sex-limited treatments and the control evolutionary treatment. c and d depict the same things as a and b, except for estimated mean male fitness.